{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'sympy'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-1-ec3cbf0d549f>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mimport\u001b[0m \u001b[0msympy\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0msp\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mnumpy\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'sympy'"
     ]
    }
   ],
   "source": [
    "import sympy as sp\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [],
   "source": [
    "# basis vectors in x, y, and t\n",
    "x=sp.Symbol('x')\n",
    "y=sp.Symbol('y')\n",
    "t=sp.Symbol('t')\n",
    "#model and data amplitudes\n",
    "Ad=sp.Symbol('A_d')\n",
    "Am=sp.Symbol('A_m')\n",
    "\n",
    "# domain width\n",
    "L=sp.Symbol('L', positive=True, finite=True)\n",
    "# domain time dinension\n",
    "T=sp.Symbol('T', positive=True, finite=True)\n",
    "\n",
    "# data errors:\n",
    "s=sp.Symbol('\\sigma')\n",
    "# data density:\n",
    "rd=sp.Symbol('rho', positive=True, finite=True)\n",
    "# wavenumbers in the x, y, and t dimensions\n",
    "kx=sp.Symbol('k_x', integer=True, positive=True)\n",
    "ky=sp.Symbol('k_y', integer=True, positive=True)\n",
    "kt=sp.Symbol('k_t', integer=True, positive=True)\n",
    "# Weighting on d3z/dx2dt\n",
    "Wxx=sp.Symbol('W_{xx}')\n",
    "# Weighting on d2x/dt2\n",
    "Wtt=sp.Symbol('W_{tt}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {},
   "outputs": [],
   "source": [
    "# basis functions:\n",
    "bf =sp.sin(2*kx*sp.pi*x/L) * sp.sin(2*ky*sp.pi*y/L) * sp.sin(2*kt*sp.pi*t/T) \n",
    "# model is the basis function times Am\n",
    "m=Am * bf\n",
    "# data is basis function times Ad\n",
    "d=Ad * bf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 195,
   "metadata": {},
   "outputs": [],
   "source": [
    "# derivatives of the basis functions\n",
    "D2xtm = sp.diff(sp.diff(sp.diff(m, x),x),t)\n",
    "D2ytm = sp.diff(sp.diff(sp.diff(m, y),y),t)\n",
    "Dxytm = sp.diff(sp.diff(sp.diff(m, x),y),t)\n",
    "D2tm = sp.diff(sp.diff(m, t),t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Combined residual:\n",
    "R=sp.integrate(\\\n",
    "               sp.integrate(\\\n",
    "    sp.integrate( rd*((d-m)/s)**2 + (Wxx*D2xtm)**2 + (Wxx*D2ytm)**2 + 2*(Wxx*Dxytm)**2+(Wtt*D2tm)**2, (x, 0, L)),\\\n",
    "    (y, 0, L)), \n",
    "               (t, 0, T))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 197,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{A_{d} L^{4} T^{4} \\rho}{L^{4} T^{4} \\rho + 16 \\pi^{4} L^{4} W_{tt}^{2} \\sigma^{2} k_{t}^{4} + 64 \\pi^{6} T^{2} W_{xx}^{2} \\sigma^{2} k_{t}^{2} k_{x}^{4} + 128 \\pi^{6} T^{2} W_{xx}^{2} \\sigma^{2} k_{t}^{2} k_{x}^{2} k_{y}^{2} + 64 \\pi^{6} T^{2} W_{xx}^{2} \\sigma^{2} k_{t}^{2} k_{y}^{4}}$"
      ],
      "text/plain": [
       "A_d*L**4*T**4*rho/(L**4*T**4*rho + 16*pi**4*L**4*W_{tt}**2*\\sigma**2*k_t**4 + 64*pi**6*T**2*W_{xx}**2*\\sigma**2*k_t**2*k_x**4 + 128*pi**6*T**2*W_{xx}**2*\\sigma**2*k_t**2*k_x**2*k_y**2 + 64*pi**6*T**2*W_{xx}**2*\\sigma**2*k_t**2*k_y**4)"
      ]
     },
     "execution_count": 197,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# solve for model amplitude that minimizes the combined residual:\n",
    "A_best=sp.solve(sp.diff(R, Am), Am)[0]\n",
    "A_best"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 1 + \\frac{16 \\pi^{4} W_{tt}^{2} \\sigma^{2} k_{t}^{4}}{T^{4} \\rho} + \\frac{64 \\pi^{6} W_{xx}^{2} \\sigma^{2} k_{t}^{2} k_{x}^{4}}{L^{4} T^{2} \\rho} + \\frac{128 \\pi^{6} W_{xx}^{2} \\sigma^{2} k_{t}^{2} k_{x}^{2} k_{y}^{2}}{L^{4} T^{2} \\rho} + \\frac{64 \\pi^{6} W_{xx}^{2} \\sigma^{2} k_{t}^{2} k_{y}^{4}}{L^{4} T^{2} \\rho}$"
      ],
      "text/plain": [
       "1 + 16*pi**4*W_{tt}**2*\\sigma**2*k_t**4/(T**4*rho) + 64*pi**6*W_{xx}**2*\\sigma**2*k_t**2*k_x**4/(L**4*T**2*rho) + 128*pi**6*W_{xx}**2*\\sigma**2*k_t**2*k_x**2*k_y**2/(L**4*T**2*rho) + 64*pi**6*W_{xx}**2*\\sigma**2*k_t**2*k_y**4/(L**4*T**2*rho)"
      ]
     },
     "execution_count": 198,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# simplify the denominator:\n",
    "sp.expand(1/(A_best/(Ad)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This expression is equivalent to A_best"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{A_{d} L^{4} T^{4} \\rho}{L^{4} T^{4} \\rho + 16 \\pi^{4} \\sigma^{2} k_{t}^{2} \\left(L^{4} W_{tt}^{2} k_{t}^{2} + 4 \\pi^{2} T^{2} W_{xx}^{2} \\left(k_{x}^{2} + k_{y}^{2}\\right)^{2}\\right)}$"
      ],
      "text/plain": [
       "A_d*L**4*T**4*rho/(L**4*T**4*rho + 16*pi**4*\\sigma**2*k_t**2*(L**4*W_{tt}**2*k_t**2 + 4*pi**2*T**2*W_{xx}**2*(k_x**2 + k_y**2)**2))"
      ]
     },
     "execution_count": 200,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.simplify(Ad*rd/s**2 /  (rd/s**2 + 16*sp.pi**4*(Wtt**2*kt**4/T**4 + 4*sp.pi**2*kt**2*Wxx**2*(kx**2+ky**2)**2/(L**4*T**2))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to look at the relationship between data bias errors and large-wavelength second derivatives, we look for the period where the second term in the denominator becomes equal to the first.  \n",
    "\n",
    "To suppress oscillations of period  < 0.125 yr, for data errors of 0.1 m, with 1 data point /3 km:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 213,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6919.4303540850005\n"
     ]
    }
   ],
   "source": [
    "tau=0.125  # 1/8 year period\n",
    "Ti=1 # 1 year\n",
    "kti = Ti/tau\n",
    "rho_tr=4/3000   # 4 measurements/year/3 km\n",
    "si=0.1 # 0.1 m error\n",
    "Wtti = np.sqrt((Ti**4*rho_tr)/(16*np.pi**4*kti**4*si**2))\n",
    "print(1/Wtti)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 216,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "26798.838526613294\n"
     ]
    }
   ],
   "source": [
    "tau=0.125\n",
    "Ti=1\n",
    "kti = Ti/tau\n",
    "si=0.1\n",
    "rho_dp=4*((2*3000/30)/3000/3000)   # 4 measurements/year*(*2 tracks*3km /60 m/pt)/(3 km x 3km)\n",
    "Wtti = np.sqrt((Ti**4*rho_dp)/(16*np.pi**4*kti**4*si**2))\n",
    "print(1/Wtti)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 225,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0003968803415078377\n"
     ]
    }
   ],
   "source": [
    "# if we want to suppress 1-km scale variations in dh/dt at a temporal scale of 0.25 yr:\n",
    "Li=1000\n",
    "wvl=500\n",
    "Ti=1\n",
    "tau=0.25\n",
    "kti=Ti/tau\n",
    "kxi=Li/wvl\n",
    "Wxxi = Li**2*Ti/(8*np.pi**3*si*kti*kxi**2)\n",
    "print(1/Wxxi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 226,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "kjer\n"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
